KrigingInitialize Subroutine

SUBROUTINE KrigingInitialize &
!
(anisotropy, nlags, maxlag, count, neighbors, nclosest)
    
IMPLICIT NONE

!Arguments with intent(in):
INTEGER (KIND = short), INTENT(IN) :: anisotropy !! if = 1, anisotropy is considered
INTEGER (KIND = short), INTENT(IN) :: nlags !! the number of discrete distances used for comparing data, if 0 it is set to default value
REAL (KIND = float), INTENT(IN) :: maxlag !! maximum distance to consider for semivariogram assessment (m), if 0 it is automatically computed
INTEGER (KIND = short), INTENT(IN) :: count !!total number of available observations
INTEGER (KIND = short), INTENT(IN) :: neighbors !! number of neighbors to consider for interpolation

!Arguments with intent(out)::
INTEGER (KIND = short), INTENT(OUT) :: nclosest

!local declarations:
INTEGER (KIND = short) :: lagnumbers !!number of lags used in computing semivariogram, default = 15
REAL (KIND = float) :: varcoverage !! variogram coverage, distance to consider for semivariogram assessment, (m), anisotropy case
REAL (KIND = float) :: varcoverageEW !! variogram coverage in East-South direction (0°)
REAL (KIND = float) :: varcoverageNESW !! variogram coverage in NorthEast-SouthWest direction (45°)
REAL (KIND = float) :: varcoverageNS !! variogram coverage in North-Sout direction (90°)
REAL (KIND = float) :: varcoverageNWSE !! variogram coverage in NorthWest-SoutEast direction (315°)
REAL (KIND = float) :: width (4)  !! distance between lags (m)
INTEGER (KIND = short) :: i, j, m, k

!----------------------------------end of declarations-------------------------

!set number of lags and allocate vectors
IF (nlags > 0) THEN
    lagnumbers = nlags
ELSE
    lagnumbers = 15 !default
END IF


IF (anisotropy == 1) THEN 
    ndir = 4
ELSE !only one direction is considered
    ndir = 1   
END IF

!allocate memory for semivariogram
ALLOCATE ( dir (ndir) ) !angle between the direction of the variogram and the x axis
ALLOCATE ( lagNumber (ndir) ) !number of lags for each direction
ALLOCATE ( dist (ndir, lagnumbers) ) !distance for each direction and lag
ALLOCATE ( empVariogram (ndir, lagnumbers) ) !experimental variogram for each direction and lag
ALLOCATE ( modVariogram (ndir, lagnumbers) ) !modelled variogram for each direction and lag
ALLOCATE ( pairNumber (ndir, lagnumbers) ) !number of pairs for each direction and lag

!set distance for semivariogram computation
IF (maxlag > 0.) THEN
    varcoverage = maxlag
ELSE !set to maximum distance between points
    IF (anisotropy == 0) THEN
         varcoverage = 2. ** 0.5 * MAXVAL (dist_pairs) / 2.  !sqrt(2) *1/2 Maximum point distance
    ELSE !search for maximum distance along  directions
       varcoverageEW   = 0.
       varcoverageNESW = 0.
       varcoverageNS   = 0.
       varcoverageNWSE = 0.
       DO j = 1, count !loop through pairs
           DO k = 1, count
               !EW direction
               IF ( dir_pairs (j,k) < pi / 8. .OR. dir_pairs (j,k) >= 15./8. * pi ) THEN 
                   IF ( dist_pairs (j,k) > varcoverageEW ) varcoverageEW = dist_pairs (j,k)
               END IF
               
               !NE-SW direction
               IF ( dir_pairs (j,k) >= pi / 8. .AND. dir_pairs (j,k) < 3./8. * pi  ) THEN 
                   IF ( dist_pairs (j,k) > varcoverageNESW ) varcoverageNESW = dist_pairs (j,k)
               END IF
               
               !N-S direction
               IF ( dir_pairs (j,k) >= 3./8. * pi .AND. dir_pairs (j,k) <= 1./2. * pi .OR. &
                    dir_pairs (j,k) > 3./2. * pi .AND. dir_pairs (j,k) < 13./8. * pi  ) THEN 
                   IF ( dist_pairs (j,k) > varcoverageNS ) varcoverageNS = dist_pairs (j,k)
               END IF
                    
               !NW-SE direction
               IF ( dir_pairs (j,k) >= 13./8. * pi .AND. dir_pairs (j,k) < 15./8. * pi  ) THEN 
                   IF ( dist_pairs (j,k) > varcoverageNWSE ) varcoverageNWSE = dist_pairs (j,k)
               END IF     
               
           END DO
       END DO
                          
    END IF
END IF

varcoverageEW = 2. ** 0.5 * varcoverageEW / 2.
varcoverageNESW = 2. ** 0.5 * varcoverageNESW / 2.
varcoverageNS = 2. ** 0.5 * varcoverageNS / 2.
varcoverageNWSE = 2. ** 0.5 * varcoverageNWSE / 2.



!set distance between lags (m). Assume even spaced lags
IF (anisotropy == 0) THEN
    width = varcoverage / lagnumbers
ELSE
    width (1) = varcoverageEW / lagnumbers
    width (2) = varcoverageNESW / lagnumbers
    width (3) = varcoverageNS / lagnumbers
    width (4) = varcoverageNWSE / lagnumbers
END IF

!populate lags vectors
DO i = 1, ndir
    DO j = 1, lagnumbers
       dist (i,j) = width (i) * j
    END DO
END DO

!set tolerance
IF (anisotropy == 0) THEN
  lagTolerance = width (1) / 2. !by default tolerance is half lag width
ELSE
    DO i = 1, ndir
        lagTolerance (i) = width (i) / 2.
    END DO
END IF
!initialize
DO i = 1, ndir !number of lags is the same for each direction
  lagNumber (i) = lagnumbers
END DO

IF (anisotropy) THEN
    !set 4 directions: E-W (0°), NE-SW (45°), N-S (90°), NW-SE (315°)
    dir (1) = 0. !E-W
    dir (2) = 45. * degToRad ! NE-SW
    dir (3) = 90. * degToRad ! N-S
    dir (4) = 315. * degToRad ! NW-SE
ELSE
    !direction is not used
    dir (1) = 0. 
END IF

pairNumber = 0
empVariogram = 0.

ieva = 0
nugget = 0.0
range = 0.0
partialSill = 0.0
anisotropyAngle = 0.0
anisotropyRatio = 1.0
objectiveFunction = 4 !set objective function to minimize

npep = 1 !1 = nugget is estimated
nvar = 0 !1 means variance = experimental variance

!find the dimension of matrixes for interpolation
IF (neighbors > 0) THEN
    IF ( neighbors > count ) THEN !too many neighbors respect to available observations
         nclosest = count !limit dimension to available observations number
    ELSE
         nclosest = neighbors
    END IF
ELSE !neighbors = 0 => include all observations
    nclosest = count
END IF

!Allocate memory for interpolation
ALLOCATE ( weights ( nclosest + 1 ) )
ALLOCATE ( hest ( nclosest ) )
ALLOCATE ( cvest ( nclosest + 1 ) )
ALLOCATE ( controlpts ( nclosest ) )
ALLOCATE ( hobs ( nclosest,nclosest ) )
ALLOCATE ( cvobs( nclosest + 1, nclosest + 1 ) )
ALLOCATE ( last (nclosest) )


RETURN
END SUBROUTINE KrigingInitialize